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© Multiple processing of radiographic images based on a pyramidal image decomposition. 

© A method of multiple processing a digital repre- 
sentation of a radiographic image is disclosed that is 
based on a pyarmidal decomposition of said image. 
A pyramidal decomposition is stored and retrieved to 
be applied to at least two different processing cy- 
cles, processed images are obtained by application 
of a reconstruction algorithm. 
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Field of the invention. 

The present invention relates to a method and 
apparatus for producing multiple differently pro- 
cessed digital images. 

More in particular the invention relates to such 
a method for use in a medical radiographic imaging 
system, such as a computed radiography system 
or a computed tomography system. 

Background of the invention 

In the field of digital radiography a wide variety 
of image acquisition techniques have been devel- 
oped rendering a digital representation of a radiog- 
raphic image. 

Among such techniques are computerised 
tomography, nuclear magnetic resonance, ultra- 
sound detection, detection of a radiation image by 
means of a CCD sensor or a video camera, radiog- 
raphic film scanning etc. 

In still another technique a radiation image, for 
example an X-ray image of an object, is stored in a 
screen comprising a photostimulable phosphor 
such as one of the phosphors described in Eu- 
ropean patent publication 503 702, published on 
16.09.92 and US Ser. 07/842,603. The technique 
for reading out the stored radiation image consists 
of scanning the screen with stimulating radiation, 
such as laser light of the appropriate wavelength, 
detecting the light emitted upon stimulation and 
converting the emitted light into an electric repre- 
sentation for example by means of a photomul- 
tipiier and finally digitizing the signal. 

One of the benefits of a digital radiographic 
system resides in the possibility of processing the 
digital image representation before display or hard 
copy recording. The term "processing" in this con- 
text means any kind of image-processing such as 
noise filtering, contrast enhancement, data com- 
pression etc. 

In some digital imaging systems the same 
original has to be processed in different ways to 
produce multiple hard-copies or displays originat- 
ing from the same image. These different versions 
of one image may be helpful for a radiologist in 
making a specific diagnosis. 

For example different versions of one original 
image may be generated that are processed taking 
into account different contrast enhancing modifying 
curves or different window level settings etc. Time- 
consuming, mostly convolution or non-linear neigh- 
bourhood operations, such as unsharp masking for 
edge enhancement, have to be repeated if e.g. 
several images processed according to different 
kernel sizes have to be produced. 

A digital radiographic image is commonly re- 
presented by about 10 MB digital data. The com- 



putation time required for processing such an 
amount of data may extend to orders of minutes 
per processing cycle. 

Hence, in case multiple differently processed 

6 versions of one image are required, it is highly 
desirable to optimize the processing procedure as 
far as the computation time is concerned. 

Recently a new image processing technique 
has been developed. According to this technique 

/o an image (more specifically a digital signal repre- 
sentation of an image) is first decomposed into a 
multiresolution representation which represents 
localised image detail at multiple scales. For exam- 
ple, the image is decomposed into a sequence of 

rs detail images at' multiple resolution levels and a 
residual image^ at a resolution lower than the mini- 
mum of said multiple resolution levels. Next, the 
pixel values of said multiresolution representation 
are modified by means of modifying function. And 

20 finally a processed image is computed by applying 
a reconstruction algorithm to the modified multiple 
resolution representation, the reconstruction algo- 
rithm being such that if it were applied to the 
unmodified multiresolution representation then said 

25 original image or a close approximation thereof 
would be obtained. 

In a preferred embodiment the decomposition 
is pyramidal, meaning that the number of pixels in 
the components at successive resolution levels de- 

30 creases. 

Objects of the invention 

It is an object of the present invention to pro- 
as vide a method of obtaining differently processed 
image versions originating from a single radiog- 
raphic original image in a fast and computationally 
inexpensive way. 

It is a further object to provide such a method 
40 in a system wherein a radiographic image is stored 
in a photostimulable phosphor screen and wherein 
a digital representation of said image is obtained 
by reading out said phosphor screen. 

Still further objects will become apparent from 
45 the description given hereinbelow. 

Statement of the invention 

To achieve the above objects the present in- 
50 vention provides a method of processing a digital 
representation of an original radiographic image, 
comprising steps of: 

1 )-transforming said image into a multiresolution 
representation which represents localised image 
55 detail at multiple scales, 

2)-storing said multiresolution representation into 
a memory, 
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3) -producing at least two differently processed 
images by applying at least two individual pro- 
cessing cycles, each processing cycle compris- 
ing the steps of; 

4) -retrieving said multiresolution representation 5 
from said memory, 

5) -modifying the multiresolution representation 
at at least one resolution level according to a 
non-identity function of a neighbourhood of re- 
trieved values, said neighbourhood consisting of w 
values of the same resolution level which cor- 
respond to a spatially coherent region of pixels 

in said image, 

6) -obtaining said processed image by applying 

the inverse of said transform to the modified is 
multiresolution representation. 

The term memory is meant to include both a 
working memory and medium- or long term mass 
storage device whereby in case of working mem- 
ory the processing is on-line processing imme- 20 
diately following image acquisition and in case of 
medium or long time mass storage the invention 
relates to future processing. 

In the statement of the invention and the de- 
scription hereinbelow interactions performed on an 25 
image or on a so-called detail image are to be 
interpreted as meaning interactions performed on 
the digital signal representation thereof. 

It is preferred that the decomposition of the 
original image into a multiresolution representation 30 
is pyramidal so that the number of pixels in detail 
images at coarser successive resolution levels de- 
creases thereby decreasing the number of pixels to 
be processed at successive resolution levels in the 
pyramid. 35 

The decomposition of an image into such a 
multiresolution pyramidal representation has been 
described extensively in our copending European 
application 91202079.9 filed on 14 August 1991. 

In one embodiment the multiresolution repre- 40 
sentation is a sequence of detail images at multiple 
resolution levels and a residual image at a resolu- 
tion lower than the minimum of said multiple reso- 
lution levels. 

Preferably the detail image at the finest resolu- 45 
tion level is obtained as the pixelwise difference 
between the original image and an image obtained 
by low pass filtering the original image, wherein 
successive coarser resolution level detail images 
are obtained by taking the pixelwise difference 50 
between two low pass filtered versions of the origi- 
nal image, the second filter having a smaller band- 
width than the former. 

in generai the decomposition is 10 be per- 
formed so that each pixel value in said original 55 
image is equal to the sum of the corresponding 
pixel value of said residual image incremented by 
the corresponding pixel value of each of said detail 



images, said residual and detail images being 
brought into register . with the original image by 
proper interpolation if their number of pixels is not 
equal to the number of pixels of the original image, 
and so that 

i) the mean of all pixel values in every detail 
image is zero; 

ii) the spatial frequency of every detail image" is 
limited to a specific frequency band, said fre- 
quency band being defined as the compact re- 
gion in the spatial frequency domain which con- 
tains nearly alt (say 90%) of the spectral energy 
of the basic frequency period of said discrete 
detail image, adjusted to the original spatial fre- 
quency scale- if said detail image contains less 
pixels than said original image; 

iii) every detail image corresponds to a different 
spatial frequency band, in such a way that the 
entire spatial frequency domain ranging from -pi 
to pi radians per pixel along both spatial fre- 
quency axes is covered by said spatial fre- 
quency bands associated with all said detail 
images considered within the decomposition; 

iv) each spatial frequency band associated with 
one of said detail images may partially overlap 
the neighbouring bands without being fully in- 
cluded by a frequency band associated with 
another detail image; 

v) the number of pixels within each detail image 
is at least the number of pixels required by the 
Nyquist sampling criterion, so as to avoid al- 
iasing, vi) at least two of said spatial frequency 
bands are considered in the course of said 
decomposition. 

Preferably the detail images at successively 
coarser resolution levels are obtained as the result 
of each of K iterations of the following steps: 

a) computing an approximation image at a next 
coarser level by applying a low pass filter to the 
approximation image corresponding to the cur- 
rent iteration, and subsampling the result in pro- 
portion to the reduction in spatial frequency ban- 
dwidth, using however the original image as 
input to said low pass filter in the course of the 
first iteration; 

b) computing a detail image as the pixelwise 
difference between the approximation image 
corresponding to the current iteration and the 
approximation image at a next coarser resolution 
level computed according the method sub 4.a), 
both images being brought into register by prop- 
er interpolation of the latter image; 

and wherein the residual image is equal to the 

: *: : n ^J,. nn< r4 U. . iUa !t/M"t(inn 

djJfjiUAimciuuii "nayo yj\ uuuuou uy mo iuot ugiwuvn, 

and wherein said processed image is computed by 
iterating K times the following procedure starting 
from the coarsest detail image and the residual 
image: 
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computing the approximation image at the current 
resolution level by pixelwise adding the detail im- 
age at the same resolution level to the approxima- 
tion image at the coarser resolution level corre- 
sponding to the previous iteration, both images 5 
being brought into register by proper interpolation 
of the latter image, using however the residual 
image instead of said coarser approximation image 
in the course of the first iteration. 

The subsampling factor is preferably 2 and ;o 
said low pass filter has an impulse response which 
approximates a two-dimensional gaussian distribu- 
tion. 

An apparently increased amount of pyramid 
values relative to the original image array is no 75 
problem since values of finer pyramid levels can 
be represented by fewer bits than values of coarser 
levels without affecting image quality. 

An alternative decomposition procedure com- 
prises decomposing said original image into a 20 
weighted sum of predetermined basic detail im- 
ages at multiple resolution levels and a residual 
basic image by applying a transform to said image, 
said transform yielding a set of detail coefficients 
each expressing the relative contribution to the 25 
original image of one of a set of basis functions 
representing said basic detail images and a resid- 
ual coefficient representing the relative contribution 
to the original image of a basis function represent- 
ing said basic residual image, whereby said basis 30 
functions are continuous and non-periodic and have 
zero mean value except for the basis function that 
represents the basic residual image, and wherein 
said transform is characterised in that there exists 
an inverse transform which returns the original im- 35 
age or a close approximation thereof when being 
applied to said transform coefficients. 

Preferably these basis functions are orthogonal. 
Still more preferable the functions are discrete 
wavelets. 40 

The invention is advantageous in that most 
image enhancement tasks can be done using a 
multiresolution representation of an image which, 
according to this invention, is stored and retrieved 
to be used for each processing task. 45 

Different processing cycles can be performed 
on a digital representation of a single original im- 
age in a fast and convenient way since the com- 
putationally most expensive step namely the mul- 
tiresolution decomposition of an image is only per- 50 
formed once and is retrieved for application of any 
kind of processing. 

Even the final resolution level can be selected 
in accordance with the destination of the processed 
image: iconic image, preview image, diagnostic im- 55 
age on screen or on film. 

The processing referred to comprises any kind 
of image processing such as 



- modification of detail contrast by modification 
of the values of detail images according to at 
least one non-linear monotonically increasing 
odd conversion function with a slope that 
gradually decreases with increasing argument 
values as has been described in our Eu- 
ropean application 91202079.9 filed 92.07.30 
and US Ser. 07/924,905. 

- noise reduction by attenuating pyramid val- 
ues taking into account the locally estimated 
image content, a method described exten- • 
.sively in our copending European application 
92201802.3 filed 92.06.19; 

- edge enhancement by increasing values of 
the finer resolution levels in the pyramid rela- 
tive to the intermediate resolution levels, 

- suppressing gradually evolving signal compo- 
nents across the image by decreasing the 
values of the coarser resolution levels relative 
to the intermediate levels, 

- or any combination of these processing oper- 
ations. 

- different amounts of processing can be ob- 
tained by varying processing parameters. 

Finally in all these cases a processed image is 
reconstructed by applying the inverse of the de- 
composition transformation. This inverse transfor- 
mation is also described extensively in the already 
mentioned European application 91202079.9 filed 
91.08.14 and US Ser. 07/924,905. 

The inverse transformation is such that when 
applied to all unmodified detail images and the 
residual image into which the original image has 
been decomposd, the original image or a close 
approximation thereof would result. 

Particular aspects of the present invention as 
well as preferred embodiments thereof will be illus- 
trated by means of the following drawings in which 
Fig. 1 is a general view of a system in which the 
method of the present invention can be applied, 
Fig. 2 is a detailed view of a system for reading 
an image stored in a photostimulable phosphor 
screen, 

Fig. 3a and 3b schematically illustrates the data 
processing performed on the read-out image 
signal, 

Fig. 4 illustrates a specific decomposition sys- 
tem, 

Fig. 5 is an example of a filter used in the 
decomposition procedure, 
Fig. 6 illustrates a specific reconstruction pro- 
cess, 

Fig. 7 illustrates a noise suppression process, 
Fig. 8 is a plot of a modifying function used for 
the purpose of contrast enhancement, 
Fig. 9 is a plot of a mapping curve. 
Fig. 1 shows a system in which the present 
invention can be applied. A radiation image of an 
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object was recorded on a photostimulable phos- 
phor screen (3) by exposing (2) said screen to x- 
rays transmitted through an object (not shown). 
The stimulable phosphor screen was conveyed in a 
cassette (3) provided with an electrically erasable 
programmable read-only memory (EEPROM). In an 
identification station (4) various kinds of data 
(name, date of birth etc) and data relating to the 
exposure and/or to the signal processing were re- 
corded onto the EEPROM. 

In a radiation read-out apparatus (1) the latent 
image stored on the photostimulable phosphor 
screen was read-out. Then the image signal was 
sent to the image processor (7). After processing 
the image signal was sent to an output device (6) 
more specifically a laser recorder. 

Figure 2 shows one embodiment of an image 
read-out unit. This figure shows a photostimulable 
phosphor screen (8) that has been exposed to an 
X-ray image of an object. 

In the radiation image readout apparatus the latent 
image stored in the photostimulable phosphor 
screen is read out by scanning the phosphor sheet 
with stimulating rays emitted by a laser (9). The 
stimulating rays are deflected according to the 
main scanning direction by means of a gal- 
vanometric deflection device (10). The secondary 
scanning motion is performed by transporting the 
phosphor sheet in the direction perpendicular to 
the scanning direction. A light collector (11) directs 
the light obtained by stimulated emission onto a 
photomultiplier (12) where it is converted into an 
electrical signal, which is next sampled by a sam- 
ple and hold circuit (13), and converted into a 12 
bit digital signal by means of an analog to digital 
converter (14) . The signal is also applied to a 
square root amplifier so that the output image 
representing signal also called 'original or raw* im- 
age (15) is a 12 bit signal which is proportional to 
the square root of applied exposure values and 
represents the pixel value in 2048x2496 pixels. 

From the output of the read-out apparatus the 
original image is sent to an image processor (nu- 
meral 7 in figure 1). 

The sequence of the different processing steps 
performed on the image signal is illustrated in 
figures 3a and 3b. 

The invention is based on a pyramidal de- 
composition of the image signal into detail images 
at multiple resolution levels and a residual image. 
The image components (i.e. detail images at mul- 
tiple resolution levels and residual image) are 
stored so that they can be retrieved to be sub- 
jected to at least two different types of processing 
or to at least two processing cycles of the same 
type but applied to a different degree. 

After processing the modified detail images 
and the residual image are recombined by applica- 



tion of a reconstruction algorithm to obtain a pro- 
cessed image representation that is finally con- 
verted into a visible image (display or hard copy) 
The display or hard-copy can be composed of 
5 a combination of differently processed images 
originating from one original image or alternatively 
more than one hard copy or display each repre- 
senting a differently processed version of one 6'rigi- 
nal image can be generated. 
w Each of the steps performed on the image 

signal is described in detail hereinbelow. 

- One embodiment of the decomposition process 
is depicted in Fig. 4. In the decomposition section 
the original image is filtered by means of a low 
75 pass filter 20, and subsampled by a factor of two, 
which is implemented by computing the resulting 
low resolution approximation image g1 only at ev- 
ery other pixel position of every alternate row. A 
detail image bO at the finest level is obtained by 
20 interpolating the low resolution approximation g1 
with doubling of the number of rows and columns, 
and pixelwise subtracting the interpolated image 
from the original image. 

The interpolation is effectuated by the inter- 
25 polator 21, which inserts a column of zero values 
every other column, and a row of zero values every 
other row respectively, and next convolves the ex- 
tended image with a low pass filter. The subtraction 
is done by the adder 22. 
30 The same process is repeated on the low reso- 

lution approximation g1 instead of the original im- 
age, yielding an approximation of still lower resolu- 
tion g2 and a detail image b1 . 

A sequence of detail images bi, i = 0. .L-1 and 
35 a residual low resolution approximation gL are ob- 
tained by iterating the above process L times." 

The finest detail image bO has the same size 
as the original image. The next coarser detail im- 
age b1 has only half as many rows and columns as 
40 the first detail image bO. At each step of the 
iteration the maximal spatial frequency of the re- 
sulting detail image is only half that of the previous 
finer detail image, and also the number of columns 
and rows is halved, in accordance with the Nyquist 
45 criterion. After the last iteration a residual image gL 
is left which can be considered to be a very low 
resolution approximation of the original image. In 
the extreme case it consists of only 1 pixel which 
represents the average value of the original image . 
50 The filter coefficients of the low pass filter of 

the preferred embodiment are presented in Fig. 5. 
They correspond approximately to the samples of 
a two dimensional gaussian distribution on a 5x5 
grid. The same filter coefficients are used for the 
55 low pass filters 20, 20\20", 20"'. at all scales. The 
same filter kernel with all coefficients multiplied by 
4 is also used within the interpolators 21, 21 ',21", 
21 M \ The factor of 4 compensates for the insertion 
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of zero pixel columns and rows. 

Next, the decomposed images are stored on a 
storage disc indicated by numeral 5 in figure 1 . 
At any time the decomposed image can be re- 
trieved for any kind of processing and processing 
to different degrees (i.e. the same processing types 
however with different processing parameters), for 
example for performing the processing cycles 
shown in figure 3a and described later. 

Then, each processing cycle applied to a de- 
composed image signal is terminated with recon- 
struction of a processed image by application of a 
reconstruction algorithm to the modified detail im- 
ages and the residual image. 

The reconstruction algorithm is identical for 
each of the processing cycles and is illustrated in 
figure 6. 

The residual image g L is first interpolated by 
interpolator 26 to twice its original size and the 
interpolated image is next pixelwise added to the 
detail image of the coarsest level b' L -i' using adder 
27. 

The resulting image is interpolated and added 
to the next finer detail image. When this process is 
iterated L times using the unmodified detail images 
b L -i bo then the original image 15 will result. 
When at the other hand the detail images are 
modified before reconstruction according to the 
findings of the present invention, then a contrast 
enhanced image 16 will result. The interpolators 
27, 27', 27", 27"' are identical to those used in the 
decomposition section. 

After reconstruction the image is subjected to a 
logarithmic conversion. 

In a great deal of the radiologic examination 
types the patient is protected against unnecessary 
exposure to x-rays by means of an x-ray opaque 
(collimation) material that is placed in the x-ray 
beam path for shielding the diagnostically irrelevant 
parts of the patient. 

However, the image data originating from the im- 
age part corresponding with the collimation ma- 
terial have an influence on the processing. Further- 
more, when reproduced unmodified, the part of the 
image corresponding with the collimation material 
may cause problems in the display, for example it 
may impair diagnosis of subtle lesions due to daz- 
zle since the unexposed image parts appear very 
bright. So, it is advantageous to exclude the data 
regarding the collimation material from further con- 
sideration during processing. 

Hence a method has been developed for deter- 
mining the signal/shadow boundary in an image so 
as to recognize the exact limits of the irradiation 
field, this method has been described in extenso in 
our copending European application entitled "Meth- 
od of recognising an irradiation field" and filed on 
the even day. According to this method many 



hypotheses (being a segmentation of an image into 
signal and shadow regions) as to the location of the 
signal/shadow boundary are built from combina- 
tions of intermediate level primitives. These inter- 
5 mediate level primitives are for example extended 
line segments. Each proposed hypothesis is sub- 
jected to a number of tests so as to detect and 
reject an incorrect hypothesis; non-rejected hyjto- 
theses are then ranked in order that a single can- 
70 didate may be chosen. 

For the purpose of reducing the computational 
effort, the exact location of the irradiation field is 
calculated by applying the above method to one 
the low resolution images resulting from the de- 
is composition processing described hereinbefore, 
namely on a low resolution image comprising 
256x312 pixel's (8 bit representation), this image is 
used as an operational tool for determining the 
irradiation field, furtheron called "the region of in- 
20 terest", this low resolution image serves as a 're- 
duced image version* and is stored on the system 
disc. 

The method described higher for delineating 
the image region of diagnostic interest results in an 

25 overlay image with the same number of elements 
as the low resolution image. The resulting overlay 
image is interpolated so as to represent 2048x2496 
pixels, a number equal to the number of pixels in 
the original image. The noninterpolated overlay im- 

30 age is stored for later use when determining the 
histogram of the region of interest in the processed 
image, as will be described furtheron. 

For hard copy recording or display the pro- 
cessed image is subjected to a signal-to-density 

35 conversion on the basis of a mapping curve defin- 
ing the relation between the individual signal values 
and the correspondingly envionsioned density val- 
ue. 

Parameters for defining the mapping curve are 
40 deduced from analysis of the histogram of the 
region of interest in the processed image. This 
region of interest is determined by selecting out of 
the pixels of the processed image only these pixels 
that belong to the image area defined by the over- 
45 lay image produced by a method described higher, 
the pixels of the region of interest in the processed 
image are then applied to a histogram calculation 
circuit. 

In a following processing step this histogram is 
50 analysed so as to determine the limits of the signal 
range relevant for display or reproduction. 

The analysis of the histogram is performed as 
described in our European application EP 
91203212.5 filed on 09.12.91 and US Ser. 
55 07/976.786/ 

The analysis of the histogram results in the defini- 
tion of a signal range to be extracted for further 
processing, this range is obtained by performing 
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the steps of 

- determining the maximum histogram frequen- 
cy. 

- selecting a value t smaller than the maximum 
histogram frequency, 

- determining (a) histogram peak(s) as a range 
of successive signal values having a cor- 
responding histogram frequency that is larger 
than t. 

- determining the most relevant histogram peak 
as the histogram peak for which the summa- 
tion of all histogram frequencies correspond- 
ing with signal values within said peak is 
maximum, 

- determining minimum and maximum signal 
values within said most relevant histogram 
peak, and 

- determining extreme values of the signal 
range to be extracted as said minimum value 
decreased with a small offset do and said 
maximum value increased by a small offset 
value di . 

For example do was equal to 0.2 log exposure units, 
di was equal to 0.1 log exposure units. 

Next, the extracted signal range is used in the 
process of defining the mapping curve as de- 
scribed in our European application EP 91203209.1 
filed 9 December 1991 and US Ser. 07/978,091. 
The mapping curve was determined as follows: 
First the minimum density value D smln and the 
maximum density value D smax envisioned in the 
hard copy were defined, D fimin was equal to fog 
density and D smax was equal to 3.0. These param- 
eters were obtained from a parameter table and are 
a function of the examination type. 
Then a canonical function defined in an orthogonal 
coordinate system between x Q , xi and y mim y max 
was retrieved from the internal memory of the 
signal processor. This function is also function of 
the examination type. 

Next two values S min and S max were determined 
that constitute a range wherein the conversion of 
signal values onto density values is determined by 
the specific shape of the canonical function. Signal 
values smaller than S min are mapped onto D sm i n , 
signal values greater than S max are mapped onto 

In this embodiment the latitude of said range was a 
fixed value L = 1.5 log exposure (corresponding 
with the latitude of a conventional x-ray film the 
radiologist is used to work with) and the position of 
S m in was determined relative to the diagnostically 
relevant signal range. S max was then calculated as 
S mln + L 

|-or determining the position of S min relative to the 
relevant signal range, the extreme values So and 
Si of the diagnostically relevant signal range were 
first determined by evaluation of^the image histo- 



gram. 

Then a small offset dSi = 0.3 log E was added to 
Si . This ensures that the density in the hard copy 
corresponding with the maximum value of the di- 
5 agnostically relevant signal range does not become 
too dark. 

The positioning of the range S max -S min relative to 
the range Si -So was performed by aligning a frac- 
' tion of the latter range with the same fraction of the 

10 former range. 

Mathematically this fraction can be expressed as 
A(Si +dSi-So-dS 0 ).Then the alignment can be for- 
mulated mathematically as S m i n = So +dSo + A- 
(S, +dSi-So-dS 0 ) - A.L and S max = S min + L ( A being 

15 an integer value greater than or equal to 0 and 
smaller than or equal to 1. 

Next a look' up table representing the mapping 

curve is composed and stored. 

Figure 9 shows a canonical curve used for deter- 

20 mining the mapping curve. 

Finally the mapping curve is applied to the 
region of interest of the processed image defined 
by application of the extrapolated overlay image so 
as to obtain the output image. 

25 Reference is again made to figure 3a illustrat- 

ing different types of processing as well as pro- 
cessing to different degrees that were applied to 
the retrieved detail images. 

The figure shows a first processing cycle 

30 wherein successively noise suppression, contrast 
enhancement, high frequency enhancement and 
low frequency suppression are performed prior to 
reconstruction of a processed image. 

Each of these processing sub-steps are de- 

35 scribed hereinbelow. 

An embodiment of a noise suppressing section 
comprising a section wherein the noise variance is 
estimated, is illustrated in figure 7. 
This type of noise suppression processing has 

40 been described in our copending European ap- 
plication 92201802.3 filed 92.06.19. 

Numeral 31 is a storage device wherein the 
detail images b t and the residual image g L resulting 
from the image decomposition are stored and cor- 

45 responds with the storage disc shown in fig. 3a. 
Each detail image is pixelwise transferred to a 
squaring unit 32, starling with the coarsest detail 
image. A moving average operator 33 then com- 
putes the local variance v at every pixel position by 

50 summing all squared pixels in an NxN neighbour- 
hood centered around the current target pixel (a 
neighbourhood of 15x15 elements proved to be 
adequate), and dividing the sum by the number of 
pixels in the neighbourhood. These local variance 

55 pixeis are temporarily stored in a storage device 34 
and transferred at the same time to a histogram 
computation circuit 35. A histogram is an array, the 
elements of which are called bins, each bin cor- 
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responding to a fixed sampling interval of the sig- 
nal range associated with the horizontal histogram 
axis. Each bin resides in a memory cell, all of them 
being initialised to zero before accepting the first 
pixel. For each entered variance value the histo- 
gram computation circuit selects the corresponding 
bin index and increments the associated bin value 
by one. 

After all pixels of a variance image at a particular 
resolution level have been used in this way, the 
histogram represents the occurence of every quan- 
tised variance value throughout the image. This 
local variance histogram is next supplied to a maxi- 
mum locator 36 which determines the variance 
value with the highest occurence v n in the histo- 
gram. This value is used as an estimate for the 
noise variance within the considered detail image. 
This estimated value is one of the parameters that 
is stored. The noise variance vn determined by the 
maximum locator is used as a parameter in the 
noise suppression function Svn (v), which is de- 
fined as: 

Svn (v) « 0 if v < = K * vn 
Svn (v) = 1 - K * vn/v otherwise 

where K is a fixed noise suppression factor which 
determines the amount of noise suppression to be 
applied; K = 0 implies no noise suppression. 
This function is computed and installed as a noise 
suppression look-up table 37 for every detail image 
within the decomposition. 

When a noise suppression look up table corre- 
sponding to a particular resolution level has been 
installed, all variance pixels corresponding with the 
same level are fetched from the storage device and 
tranformed into a sequence of attenuation coeffi- 
cients. The resulting pixels are computed by pixel- 
wise multiplying 38 these coefficients with the pix- 
els of the detail image at the same level, fetched 
from the storage device 31. 

This whole process is repeated for all detail images 
up to the finest level, to yield attenuated detail 
images. 

Next, the modified detail images are subjected 
to contrast enhancement as described extensively 
in our copending European applications 
91202079.9 filed 91.08.14 and US Ser. 07/924,905. 

This step is performed by modifying the pixels 
of the detail images (after noise suppression) to 
yield pixel values of a set of modified detail images 
according to at least one non-linear monotonically 
increasing odd modifying function with a slope that 
gradually decreases with increasing argument val- 
ues. 

in the modification section a lookup table is 
provided which converts every pixel value x of 
each detail image into an output value y according 



to the function: 

y = -m ' (-x/m) p if x < 0 
y = m * (x/m) p if x > = 0 

5 

where the power p is chosen within the interval 0 < 
p < 1, preferably within the interval 0.5 < p < 0.9. 
A comparative evaluation of a large number of 
computed radiography images of thorax and bones 

w by a team of radiologists indicated that p = 0.7 is 
the optimal value in most cases, m specifies the 
abscissa range : -m < = x < = m, e.g. m = 4095 if 
detail pixels are represented by 13 bits signed. 
A plot of the above function is presented in Fig. 8. 

75 If depending on the used decomposition method 
either the pixels of each detail image or otherwise 
the detail weighting coefficients as obtained from 
one of the above decomposition methods, are con- 
verted according to the above function, then all 

20 details with a low amplitude will be boosted relative 
to the image details wich originally have a good 
contrast. In this respect the above power function 
proved to perform very well, but it is clear that an 
infinite variety of monotonically increasing odd 

25 mapping functions can be found that will enhance 
subtle details. The main requirement is that the 
slope of said mapping function is steeper in the 
region of argument values that correspond to small 
detail image pixel values or coefficient values than 

30 it is in the region of large detail pixel or coefficient 
values. 

In an alternative embodiment excessive noise am- 
plification can be avoided by using a composite 
mapping function: 

35 

y = -m * (-x/m) p2 if -m < = x < -c 
y = -m * ( c/m) p2 * (-x/c) p1 if -c < = x < 0 
y = m * ( c/m) p2 * ( x/c) pl if 0 < = x < c 
y = m * ( x/m) p2 ifc< = x< = m 

40 

where the power p2 is chosen within the interval 0 
< p2 < 1, preferably 0.5 < p2 < 0.9, and most 
preferably p2 = 0.7 (however the preferred value 
of p2 depends upon the kind of radiological exami- 

45 nation), where the power p1 is not smaller than p2: 
p1 > = p2, where the cross-over abscissa c speci- 
fies the transition point between both power func- 
tions: 0 < c < m, and preferably c is very small 
relative to m; and where m specifies the abscissa 

50 range: 

-m < = x < = m. 

Then after the contrast enhancement step the 
55 modified detail images were subjected to a high 
frequency emphasising step. Edge enhancement is 
obtained by increasing values of the finer resolution 
levels in the pyramid relative to the intermediate 
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resolution levels. 

The following processing step comprises sup- 
pressing gradually evolving signal components 
across the image by decreasing the values of the 
coarser resolution levels relative to the intermediate 
levels. 

Finally the processed detail images and the 
residual image were subjected to a reconstruction 
algorithm as described higher and to the recon- 
structed processed image a logarithmic conversion 
was applied. 

Then the reconstructed image was converted 
into a visible image. 

Apart from the above described processing cy- 
cle, figure 3a shows alternative processing cycles 
that can be performed on the same detail images 
into which the original image has been decom- 
posed and that are stored. 

In this example the other cycles are either 
alternative combinations of the above described 
processing steps for example only comprising con- 
trast enhancement and no noise suppression nor 
high frequency emphasis nor low frequency sup- 
pression. 

Still other cycles illustrate identical combina- 
tions of processing steps, however being applied to 
a different degree. For example it is possible to 
amend parameters of the noise suppression func- 
tion Svn (y), which is defined as: 

Svn (v) = 0 if v < = K " vn 
Svn (v) = 1 - K * vn/v otherwise 

where K is a fixed noise suppression factor which 
determines the amount of noise suppression to be 
applied. 

Or, alternatively to amend parameters in the 
modifying function applied for contrast enhance- 
ment which is given by 

y = -m * (-x/m) p if x < 0 
y » m * (x/m)P if x > = 0 

where the power p is chosen within the interval 0 < 
P< 1.. 

It will be clear that many alternatives may be 
envisioned that fall within the scope of the present 
invention. The invention is not limited to the de- 
scribed processing steps, combinations or process- 
ing degrees. 

Claims 

1- A method of processing a digital representa- 
tion of an original radiographic image, compris- 
ing steps of: 

1)-transforming said image into a multi- 
resolution representation which represents 



localised jmage detail at multiple scales, 

2) -storing said multiresolution representation 
into a memory, 

3) -producing at least two differently pro- 
5 cessed images by applying at least two 

individual processing cycles, each process- 
ing cycle comprising the steps of: 

4) -retrieving said multiresolution representa- 
tion from said memory, 

io 5)-modifying the multiresolution representa- 

tion at at least one resolution level accord- 
ing to a non-identity function of a neigh- 
bourhood of retrieved values, said neigh- 
bourhood consisting of values of the same 

75 resolution level which correspond to a spa- 

tially coherent region of pixels in said im- 
age, 

6)-obtaining a processed image by applying 
the inverse of said transform to the modified 
20 multiresolution representation. 

2. A method according to claim 1 wherein the 
multiresolution representation is a pyramidal 
representation wherein the number of pixels in 

25 each detail image decreases at each coarser 

resolution level. 

3. A method according to claim 1 wherein said 
multiresolution representation is obtained as a 

30 sequence of detail images at multiple resolu- 

tion levels and a residual image at a resolution 
level lower than the minimum of said multiple 
resolution levels. 

35 4. A method according to claim 3 wherein the 
detail image at the finest resolution level is 
obtained as the pixel wise difference between 
the original image and an image obtained by 
low pass filtering the original image, and 

40 wherein the successive coarser resolution level 

detail images are obtained by taking the pixel- 
wise difference between two low pass filtered 
versions of the original image, the second hav- 
ing a smaller bandwidth than the former. 

45 

5. A method according to claim 3 

- wherein the detail images at successive 
coarser resolution levels are obtained as 
a result of K iterations of the following 

so steps: 

a) computing an approximation image at 
a next coarser level by applying a low 
pass filter to the approximation image 
corresponding to the current iteration, 

55 ans subsampling the result in proportion 

to the reduction in spatial frequency ban- 
dwidth, using the original image as input 
to said low pass filter in the course of the 



9 



17 



EP 0 610 604 A1 



18 



first iteration; 

b) computing a detail image as the pixel- 
wise difference between the approxima- 
tion image corresponding to the current 
iteration and the approximation image at 
a next coarser resolution level computed 
according to the method sub 4a;both im- 
ages being brought into register by prop- 
er interpolation of the latter image; 

- and wherein the residual image is equal 
to the approximation image produced by 
the last iteration, 

- and wherein said processed image is 
computed by iterating K times the follow- 
ing procedure starting from the coarsest 
detail image and the residual image: 
computing the approximation image at 
the current resolution level by pixelwise 
adding the modified detail image at the 
same resolution level to the approxima- 
tion image at the coarser resolution level 
corresponding to the previous iteration, 
both images being brought into register 
by proper interpolation of the latter im- 
age, using the residual image instead of 
the coarser approximation image in the 
course of the first iteration. 

6. A method according to claim 5 wherein said 
subsampling factor is 2, and said low-pass 
filter has an impulse response which approxi- 
mates a two-dimensional gaussian distribution. 

7. A method according to claim 1 wherein said 
multiresolution representation comprises detail 
transform images at multiple resolution levels 
and a residual coefficient, each detail image 
comprising a set of transform coefficients ex- 
pressing the relative contribution to the original 
image of the corresponding basis function out 
of a set of predetermined basis functions, each 
of said functions representing local detail at a 
specific resolution level and being non-peri- 
odic, and having zero mean value, and wherein 
said transform is characterised by the exis- 
tence of an inverse transform which returns the 
original image or a close approximation thereof 
when being applied to said transform coeffi- 
cients and said residual coefficient. 

8- A method according to claim 7 wherein said 
basis functions are orthogonal. 

9= A method according to claim 7 wherein said 
basis functions are discrete wavelets. 



ble phosphor screen and wherein said digital 
image representation is obtained by scanning 
said screen with stimulating irradiation, detect- 
ing light emitted upon stimulation and convert- 
5 ing said detected light into a digital representa- 

tion. 

11. A method according to claim 1 wherein satd 
function depends on the value of a corre- 

jo sponding pixel in said original image. 

12. A method according to claim 1 wherein said 
function has a slope that gradually decreases 
with increasing argument values with the ex- 

75 ception of a region of lowest absolute argu- 

ment values where the slope is constant or 
increasing' 

13. A method according to claim 1 wherein said 
20 function is not identically defined at each reso- 
lution level in said decomposition, such that 
modifications applied to finer resolution levels 
ampplify image detail more than modifications 
applied to intermediate resolution levels. 

25 

14. A method according to claim 1 wherein said 
function is not identically defined at each reso- 
lution level in said decomposition, such that 
modifications applied to coarser resolution lev- 

30 els amplify image detail less than modifications 

applied to intermediate resolution levels. 

15. A method according to claim 1 wherein said 
values are modified as a function of the es- 

35 timated amount of image content and in accor- 

dance with an estimated noise level. 

16. A method according to claim 15 wherein said 
noise level is determined as the estimated 

40 noise level in each detail image. 



45 
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17. 



A method according to claim 15 wherein a set 
of noise suppression functions are computed , 
each being associated with one of said detail 
images, said functions being monotonically 
non-decreasing in one variable and para- 
metrically depending on said noise variance in 
a non-increasing monotonic way and said sup- 
pression functions being positive and asymp- 
totically reaching a maximum value equal to 
one and wherein each detail image is attenu- 
ated by multiplying it with an associated noise 
suppression function evaluated at an abscissa 
equal to said local detail image variance. 



10. A method according to claim 1 wherein said 
radiographic image is stored in a photostimula- 
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